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ABSTRACT 

We study the evolution of binary stars in globular clusters using a novel approach combining a state-of-the-art pop- 
ulation synthesis code with a simple treatment of dynamical interactions in the dense cluster core. We find that 
the combination of stellar evolution and dynamical interactions (binary-single and binary-binary) leads to a rapid 
(•f) depletion of the binary population in the cluster core. The maximum binary fraction today in the core of a typical 
f^*) ' dense cluster like 47 Tuc, assuming an initial binary fraction of 100%, is only about 5%. We show that this is in good 
\ agreement with recent HST observations of close binaries in the core of 47 Tuc, provided that a realistic distribution 
CNI . of binary periods is used to interpret the results. Our findings also have important consequences for the dynamical 
' modeling of globular clusters, suggesting that "realistic models" should incorporate much larger initial binary fractions 
D \ than has usually been done in the past. 

f"*^ ■ Subject headings: binaries: close — binaries: general — methods: n-body simulations — globular 

\ clusters: general — globular cluster: individual (NGC 104, 47 Tucanae) — stellar 

dynamics 
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1. INTRODUCTION 

Binary stars play a fundamental role in the dynam- 
ical evolution of globular clusters, providing (through 
inelastic collisions) the source of ener gy that sup- 
ports them against grayothermal collapse llFregeau et alJ 
I2003bt IC^o~et^ll99lT: ICoodman fc Hud 119891) . In the 
"binary burning" phase, a cluster can remain in quasi- 
thermal equilibrium with nearly constant core density 
and velocity dispersion for many relaxation times, in 
much the same way that a star can maintain itself in 
thermal equilibrium for many Kelvin-Hclmholtz times by 
burning hydrogen in its core. 

At present, there are very few direct measurements of 
binary fractions in clusters. However, even early obser- 
vations showed that binary fractions in globular cluster 
cores are smaller than in the solar neighborhood (e.g., 
iCote et ail l)1996|) L Recent Hubble Space Telescope ob- 
servations have provided further constr aints on the bi- 
nary fractions in many globular c lusters l)Bellazzini et all 
120021 iRubenstein fc Bailvn1ll997[) . The measured binary 
fractions in dense cluster cores are found to be very 
small. As an example, the upper limit on the core bi- 
nary fraction of NGC 6397 is only 5-7% l|Cool fc Boltonl 
I2002T) . 

All dynamical interactions in dense cluster cores tend 
to destroy binaries (with the possible exception of tidal 
captures, which may form binaries, but turn out to play 
a negligible role; see § 3) . Soft binaries can be disrupted 
by any tidal interaction with another passing star or bi- 
nary. Even hard binaries can be destroyed in resonant 
binary-binary encounters, which typically eject two sin- 
gle st ars and leave only one binary remaining ((Mikkolal 
Il983[). or produce physical stellar collisions and mergers 
ijBacon et al~lll996t iFreeeau et alJl2003a|) . 

In addition, many binary stellar evolution processes 
can lead to disruptions (e.g., following a supernova explo- 
sion of one of the stars) or mergers (e.g., following a com- 
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mon envelope phase). These evolutionary destruction 
processes can also be enhanced by dynamics (e.g., more 
common envelope systems form as a result of exchange 
interactions; see Rasio, Pfahl, & Rappaport 2000). 

It is therefore natural to ask whether the small binary 
fractions measured in old globular clusters today result 
from these many destruction processes, and what the ini- 
tial binary fraction must have been to explain the current 
numbers. We address these questions in this paper, by 
performing calculations that combine binary star evolu- 
tion with a treatment of dynamical interactions in dense 
cluster cores. 

2. METHODS AND ASSUMPTIONS 

Our initial conditions are described by the following 
parameters: total number of stars (single or in a bi- 
nary), N, initial mass function (IMF), binary fraction, 
fb, distribution of binary parameters (period, P, eccen- 
tricity, e, and mass ratio, q). We typically adopt stan- 
dard choices used in population synthesis studies, which 
are based on available observations for stars in the field 
and in young star clusters (Sills et al. 2003). For the 
calculations reported here, we u se the fo l lowing initial 
conditions. Wc adopt the IMF of lKroupal (|2002f) . which 
can be written as a broken power law dN oc m~ a dm, 
where a = 0.3 for 0.01 < m/M Q < 0.08, a = 1.3 for 
0.08 < tu/Mq < 0.50, a = 2.3 for m/M Q > 0.5 Wc 
consider the mass range 0.05 Mq to 100 M©. The initial 
average stellar mass is then (m) = 0.48 Mq. The binary 
mass ratio, q, is assumed to be distributed uniformly in 
the range < q < 1. This is in agre ement with observa- 
tions for q > 0.2 ijWoitas et al.l200iT) . The binary period, 
P, is taken from a uniform distribution in log 10 P over 
the range P = 0.1-10 7 d. The binary eccentricity, e, fol- 
lows a thermal distribution truncated such that there is 
no contact binary. 

We evolve all stars (single and bi nary) using the 
popu lation synthesis code Star Track (|Belczvnski et aD 
l2002|) . The evolution of single stars i s base d on the 
analytic fits provided bv iHurlev et al. ( 2000), but in- 
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eludes a more realistic determin ation of compact object 
masses ijFrver fc Kalogeral200H) . All our calculations use 
metallicity Z = 0.001, appropriate for a cluster such 
as 47 Tuc. We treat the evolution of stellar collision 
and binary merger products following the prescription of 
iHurlev et alJ|2002T ). To evolve the cluster population, we 
consider two timesteps. One is associated with the evolu- 
tionary changes in the stellar population, Ai GV , and the 
other with the rate of encounters, Ai co n (see §3). At cv 
is defined so that no more than 2% of all stars change 
their properties (mass and radius) by more than 5%. The 
global timcstcp for the cluster evolution is taken to be 
At = min[i cv , £ coU ] 

Our modeling of the cluster dynamics is highly simpli- 
fied. We assume that the core number density, n c , and 
one-dimensional velocity dispersion, a, remain strictly 
constant throughout the evolution. These quantities are 
input parameters used to calculate dynamical interaction 
rates in the cluster core (see below). While all globular 
clusters have a ~ lOkms -1 , the core density can vary by 
several orders of magnitude. Here we set n c = 10 5 pc -3 
for most calculations, representative of a fairly dense 
cluster like 47 Tuc. In general, n c is the main "knob" 
that we can turn to increase or decrease the importance 
of dynamics. Setting n c = corresponds to a traditional 
population synthesis simulation, where all binaries and 
single stars evolve in isolation after a single initial burst 
of star formation. To model a specific cluster, we match 
its observed core luminosity volume density pi, central 
velocity dispersion, and half-mass relaxation time. 

The escape speed from the cluster core can be 
estim ated from observations as v e = 2.5 a? l)Webbinkl 
Il985|) . where 03 is the three-dimensional core velocity 
dispersion. Following an interaction or a supernova 
explosion, any object that has acquired a recoil speed 
exceeding v e is removed from the simulation. For 
computing interactions in the core, the velocities of 
all objects arc assumed t o be distrib uted according 
to a lowered Maxwellian iKind I1965J) . with f(v) = 
v 2 /a(m) 2 (exp(—1.5v 2 /a(m) 2 ) — exp(—1.5v 2 /a(m) 2 )) 
with parameters o~(m) = ((m)/m) 1 ' 2 <73 (assuming 
energy cquipartition in the core) and v c . In addition, we 
use a to impose a cut-off for soft binaries entering the 
core. Any binary with maximum orbital sp eed < O.lo a 
is immediately broken into two sing le stars lIHillsll 1990ft . 

In the presence of a broad mass spectrum, the clus- 
ter core is always dominated by the most massive ob- 
jects in the cluster, which tend to concentrate there via 
mass segregation. As stars evolve, the composition of 
the core will therefore change significantly over time. 
To model mass segregation in our simulations, we as- 
sume that the probability for an object of mass m to 
enter the core after a time t s follows a Poisson distri- 
bution, p(t s ) = (l/t sc ) exp(— t s /t sc ), where the char- 
acteristic mass -segregation timesca le is given by t sc = 
10 ({m} /m) t r h (|Fregeau et al.l2002() . Here i rn is the half- 
mass relaxation time, which we assume to be constant for 
a given cluster. We fix (m) = 2M , as this value gives, 
in our model, the best fit for the ratio of core mass to 
total cluster mass in 47 Tuc. 

All objects are allowed to have dynamical interactions 
after they have entered the cluster core. We use a simple 
Monte-Carlo prescription to decide which pair of objects 



Table 1. Reference models. 



Model logn c logt rh / b>c /0.5 /wd fb 



1 


5.0 


9.0 


0.08 


0.07 


0.04 


0.65 


D3 


3.0 


9.0 


0.21 


0.25 


0.10 


0.70 


Di 


1.0 


9.0 


0.14 


0.15 


0.07 


0.66 


D6 


6.0 


9.0 


0.03 


0.02 


0.02 


0.66 


T8 


5.0 


8.0 


0.10 


0.04 


0.03 


0.22 


T10 


5.0 


10. 


0.06 


0.06 


0.04 


0.85 


B05 


5.0 


9.0 


0.05 


0.05 


0.05 


0.34 


ND 


5.0 


0.0 


0.90 


0.77 


0.50 


0.90 



Note. — n c is the core number density in pc~ 3 (assumed fixed), 
i r h is the half-mass relaxation time in yr, /b ;C is the binary fraction 
in the core, /q,s is the binary fraction for non-degenerate stars more 
massive than 0.5 Mq, / w a is the binary fraction among white dwarfs 
and /t, is the overall binary fraction in the cluster. ND is the model 
with no dynamical interactions (field population), where all stars are 
assummcd to be in the core from the beginning. Values for all binary 
fractions arc given at 14 Gyr. 



actually have an interaction during each timestep. We 
considcrer separately binary-binary and binary-single 
interactions, as well as single-single encounters (tidal 
captures and collisions). Tid al captures are treated us- 
ing th e approach described in lPortegies Zwart fc Meinenl 
l|1993ft . If the pericenter distance is less than twice 
the sum of the stellar radii, the encounter is treated as 
a physical collision and assumed to lead to a merger. 
Each dynamical interaction involving a binary is calcu- 
lated using Fewbody, a new numerical toolkit for sim- 
ulating small-A gravitational dynamics that is partic- 
ularly suitejij£_rjerfoj™ 

tions l|Fregeau et al.l2003aHFregeau fc RappaportJ2003ft . 

A more detailed desc ription of our Mont e-Carlo proce- 
dure will be given in ijlvanova et "aT1 l200l . 

Initial conditions for all our reference models are given 
in Table 1. All models have central velocity dispersion 
a i = 10km s , and initial primordial binary fraction 
fb,o = lj except for Model B05, which has = 0.5. 
Our assumed period distribution implies that about 60% 
of primordial binaries are hard. The initial number of 
stars is N = 2.5 x 10 5 , with the cluster core containing 
about 1% of the stars initially. 

3. RESULTS 

First consider our results for a typical dense cluster, 
in Model 1. Starting with 100% binaries initially, the 
final core binary fraction (at 14 Gyr), /b, C! is only 8%. 
This is strikingly low, given that the cluster started with 
all binaries. Decreasing the initial binary fraction, /^o, 
to a more reasonable (but still large) 50% reduces /b c 
further to 5%, as shown in Model B05. The dependence 
°f fb,c on fb,o is not linear. This is mainly due to mass 
segregation: decreasing /^o also increases the ratio of 
mean binary mass to mean stellar mass in the cluster, 
thereby resulting in a higher concentration of binaries in 
the core. The majority (about 75%) of destroyed binaries 
were disrupted by close dynamical encounters (or, rarely, 
following a supernova explosion). Note that some bina- 
ries that are initially hard eventually become soft after 
undergoing significant mass loss due to stellar evolution. 
About 20% of the destroyed binaries experienced merg- 
ers, typically after significant hardening through interac- 
tions. The rate of binary destruction by mergers is about 
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No dynamics 
n r = I 5 , overall 
n„=10 3 , core 
core 
core 
core 



n c =10 4 
n c =10 5 
n =10 8 




Table 2. Models of specific clusters 
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Fig. 1. — Evolution of the core and overall binary fraction in 
Model 1. Also shown is the evolution with no dynamics (ND) and 
for different core number densities (D3, D4 and D6). 



an order of magnitude higher in this model (Model 1) 
than for the corresponding field population (Model ND). 
A few percent of the binaries lost were actually not de- 
stroyed but instead were ejected from the cluster as the 
result of strong encounters. Tidal capture did not play a 
significant role; the total number of tidal capture bina- 
ries formed during the cluster lifetime is less than 1% of 
the final number of binaries in the core. While the final 
core binary fraction is extremely low, the overall cluster 
binary fraction remains high, about 65%, even after 14 
Gyr (Fig. 1). Note, however, that the surviving binaries 
include mainly low-mass systems which never entered the 
core (about 70% of the initial binaries never entered the 
core and therefore never had a chance to interact). The 
average primary mass among binaries remaining outside 
the core at 14 Gyr is 0.2 Mq. 

Let us now compare results for different central den- 
sities, in Models 1, D3, D4 and D6. The evolution of 
/b lC for these models is shown in Figure 1. For compar- 
ison, the binary fraction for the field case (Model ND) 
is also presented. As expected, the core binary fraction 
decreases as n c increases. However, if one consider the 
binary fraction for only non-degenerate stars, its behav- 
ior is different. In particular, the final binary fraction 
for non-degenerate stars more massive than 0.5 Mq, fo.5, 
is higher than / DjC for the low-density Model D3. Thus 
/h e is decreased partially through a lower binary fraction 
of degenerate objects. Degenerate objects, compared to 
non-degenerate 0.5 — 0.9M & stars, evolved from initially 
more-massive stars and are more likely to have a more 
massive companion initially. Their binary destruction 
rate is much higher, enchanced both by stellar evolution 
(mass loss and mass transfer at more advanced evolution- 
ary stages, SN explosions in a binary), and by dynamical 
interactions (large cross-section for encounters). 

Next we examine how the half-mass relaxation time af- 
fects binary fractions (Models T8 and T10). We see that, 
surprisingly, the model with shorter relaxation time has 
a higher core binary fraction. There are two compct- 



Model 


log p° L 


logn c 


or 


log t rh 


fb,c 


/o.5 


fwd 


NGC3201 


2.73 a 


3.3 


5.2 


9.2 


0.17 


0.22 


0.12 


u) Cen 


3.37 


1.0 


16 


10 


0.13 


0.16 


0.06 


M3 


3.51 


3.9 


4.8 b 


8.5 


0.16 


0.13 


0.07 


M4 


3.82 


1.3 


4.2 


8.82 


0.12 


0.10 


0.07 


47Tuc 


4.81 


5.3 


11.5 


9.48 


0.05 


0.04 


0.03 


NGC6397 


5.68 


6.5 


4.5 


8.46 


0.03 


0.02 


0.02 



Note. — p° L is the observed core luminosity density in Lq pc ' , 
n c is the core number density in pc -3 , a\ is the velocity dispersion in 
kms — 1 , r r h is the half-mass relaxation time in yr, f^ c is the binary 
fraction in the core, /0.5 is the binary fraction for non-degenerate ob- 
jects more massive than 0.5 Mq, and / w d is the ratio of WD in binaries 
to the total number of WD. Values for all binary fractions are given at 
14 Gyr. 

a unless specified, p° L and t r h for specific clusters are taken from 
IHarrisI ft99"r1). a i fromlPrvor fc Mevlanl tl993T) . 
toubath eta!] 11991) 



47 Tuc 



47 Tuc, MS[3 

NGC 6397 

NGC 3201 




Fig. 2. — The binary period distribution in the models of 47 Tuc, 
tjCen, and NGC 6397 at 14 Gyr. The shaded histogram shows the 
period distribution of binaries containing two MS stars with masses 
greater than 0.25 Mq in 47 Tuc. N[, is the total number of binaries 
and Nb, p is the number of binaries in the given period range. 



ing mechanisms that play a role here: mass segregation, 
which brings binaries into the core, and dynamical inter- 
actions, which destroy binaries in the core. A shorter seg- 
regation time increases the rate at which binaries enter 
the core but also allows less massive binaries to interact. 
Therefore, the average mass of a binary in the core will 
be smaller. However, the average time spent by a binary 
in the core also increases, so more can be destroyed, and 
the more massive binaries tend to be destroyed first as 
they have a larger interaction cross section. As a result, 
in Model T8 /b, c is higher than in T10, although the 
binary fraction of more massive binaries is smaller. 

4. DISCUSSION AND COMPARISON WITH OBSERVATIONS 

We performed several simulations with parameters 
that attempt to match those of specific globular clus- 
ters in the Galaxy (Table 2). In all cases the initial bi- 
nary fraction is 100%, so our results for final core binary 



4 



fractions represent upper limits. As before we obtain 
very low values for / b . c - For example, in the 47 Tuc 
model, /b jC is only 5%. At first glance, thi s may seem to 
conflic t with observations. In particular, lAlbrow et alJ 
lj200l derive a binary fraction for the core of 47 Tuc of 
about 13%, from observations of eclipsing binaries with 
periods in the range P ~ 4-16 d. This estimate was 
based on an extrapolation assuming a period distribu- 
tion flat in log 10 P from about 2 d to 50 yr. In Figure 2 
we show the period distribution of core binaries in our 
simulation. Note that the period range of eclipsing bi- 
naries corresponds to the peak of the distribution, while 
for larger P it drops rapidly. In particular, for binaries 
with components more massive than 0.25 Mq, the num- 
ber of systems with periods in the range 2 d to 50 yr is 
about 7 times smaller than would be predicted by a dis- 
tribution flat in log 10 P. Using the observed number of 
eclipsing binaries and those from our simulation, the ad- 
justed core binary fraction from lAlbrow et all lj2001f) is 
about 4%, which is consistent with our results. Figure 2 
also shows the period distributions for models that rep- 
resent the clusters NGC 3201 and NGC 6397. For denser 
clusters, such as NGC 6397, the peak of the distribution 
shifts toward shorter periods, while for less dense clus- 
ters, such as NGC 3201, the distribution peaks at longer 
periods and is flatter. We performed 3 additional sim- 
ulations for 47 Tuc, with / bj0 = 0.25, 0.5 and 0.75. We 
found that with decreasing /b,o, the adjusted core binary 
fraction decreases and reaches, e.g., 2.5% for /b,o = 0.5. 

An alternative estimate of the binary fraction in the 
core of 47 Tuc is ba sed on observations of BY Dra stars 
ijAlbrow et al.ll2001(K Their estimated core binary frac- 
tion, which can be considered a lower limit, is approx- 
imately 0.8%, 18 times lower than the estimate based 
on eclipsing binaries. This estimate was based on 31 
BY Dra binaries and 5 eclipsing binaries oberved in the 
period range 4-10 d. We analyzed the core binary popu- 
lation in our model in order to identify BY Dra binaries. 
We adopted the standard definition for a BY Dra bi- 
nary: primary mass in the range 0.3-0.7 Mq (see, e.g., 
IBodp fe Fekell l|1977|) ') and period in the range 4-10 d (as 
for the observed sample in 47 Tuc). The ratio between 
the total number of binaries and the number of BY Dra 
systems is found to be 38, 42 and 45 for models with 



/b,o = 1.0,0.75 and 0.5, respcctivly. With 36 binaries 
observed in this period range, the total core binary frac- 
tion is 2.9%, 3.3% and 3.5% for f hfi = 1.0,0.75 and 0.5, 
respectively Based on these results, we estimate that the 
initial binary fraction should be at least 0.5, and more 
likely in the range /b,o — 0.75 — 1. 

Even more extreme results are obtained for NGC 6397. 
This cluster is classified observationally as "core col- 
lapsed" and therefore may not be well described by our 
simplified dynamical model. Nevertheless, it is useful 
to study how the binary fraction in the cluster would 
have evolved if the very high central density had been 
constant throughout the evolution. We find that the bi- 
nary fraction for stars more massive than 0.5 Mq is ex- 
tremely low: 2% in the core at 14 Gyr with 100% bina- 
ries initially. For this cluster, there is a firm upper limit 
of 3% on the core binary fraction for stars in the mass 
range 0.45 — 0.8 M Pl a nd for binary mass ratios q > 0.45 
l|Cool fe Boltonl2002|) . For the same ranges in our model 
we find 2%, in surprisingly good agreement with the ob- 
servations. 

In addition to their implications for the interpreta- 
tion of observed binary fractions in cluster cores, our 
results also have important consequences for the theo- 
retical modeling of globular clusters. Indeed, it is clear 
that "realistic" dynamical simulations of globular clus- 
ter evolution should include large populations of primor- 
dial binaries, with initial binary fractions in the range 
~ 50%-100% (sim ilar to what is usually assu med for 
the field; see, e.g., iDuauennov fc Mayo! l|1991f) ). This 
poses a particular challenge for direct A-body simula- 
tions, where the treatment of even relatively small num- 
bers of binaries can add enormous computational costs. 
For this reason, current direct A-body simulations of star 
clusters with large initial binary fractions typically have 
A too small to be considered representative of globular 
clusters (see, e.g., Portegies Zwart et al. 2003; Wilkin- 
son et al. 2003). However, approximate methods, such 
as Monte-Carlo, do not suffer from the same limitations, 
and routinely simulate clusters with reasonably large A 
(~ 3 x 10 5 ) and binary fraction (~ 30%), yet have not yet 
included advanced treatments of binary star evolution 
(see, e.g., iFreeeau et al.ll20C)3bl and references therein). 



REFERENCES 



Albrow, M. D., Gilliland, R. L., Brown, T. M., Edmonds, P. D., 
Guhathakurta, P., & Sarajedini, A. 2001, ApJ, 559, 1060 

Bacon, D., Sigurdsson, S., & Davics, M. B. 1996, MNRAS, 281, 
830 

Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407 
Bcllazzini, M., Fusi Pecci, F., Montcgriffo, P., Messineo, M., 

Monaco, L., & Rood, R. T. 2002, AJ, 123, 2541 
Bopp, B. W., & Fekel, F. 1977, AJ, 82, 490 

Cool, A. M., & Bolton, A. S. 2002, in ASP Conf. Ser. 263: Stellar 

Collisions, Mergers and their Consequences, 163 — h 
Cote, P., Pryor, C, McClure, R. D., Fletcher, J. M., & Hesser, 

J. E. 1996, AJ, 112, 574 
Dubath, P., Meylan, C, & Mayor, M. 1997, A&A, 324, 505 
Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 
Fregeau, J. M., Cheung, P., Portegies Zwart, S. F., & Rasio, F. A. 

2003a, MNRAS (submitted) 
Fregeau, J. M., Giirkan, M. A., Joshi, K. J., & Rasio, F. A. 2003b, 

ApJ, 593, 772 

Fregeau, J. M., Joshi, K. J., Portegies Zwart, S. F., & Rasio, F. A. 

2002, ApJ, 570, 171 
Fregeau, J. M., & Rappaport, S. A. 2003, in preparation 
Fryer, C. L., & Kalogera, V. 2001, ApJ, 554, 548 



Gao, B., Goodman, J., Cohn, H., & Murphy, B. 1991, ApJ, 370, 
567 

Goodman, J., & Hut, P. 1989, Nature, 339, 40 
Harris, W. E. 1996, AJ, 112, 1487 
Hills, J. G. 1990, AJ, 99, 979 

Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 
Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 
Ivanova, N. S., Belczynski, K., Fregeau, J. M., & Rasio, F. A. 2004, 

ApJ ( in preparation ) 
King, I. R. 1965, AJ, 70, 376 
Kroupa, P. 2002, Science, 295, 82 
Mikkola, S. 1983, MNRAS, 203, 1107 

Portegies Zwart, S., Hut, P., McMillan, S., & Makino, J. 2003, 

submitted to MNRAS (astro-ph/0301041) 
Portegies Zwart, S. F., & Meinen, A. T. 1993, A&A, 280, 174 
Pryor, C, & Meylan, G. 1993, in ASP Conf. Ser. 50: Structure and 

Dynamics of Globular Clusters, 357 — h 
Rubenstein, E. P., & Bailyn, C. D. 1997, ApJ, 474, 701 
Sills, A., & et al. 2003, New Astronomy, 8, 605 
Webbink, R. F. 1985, in IAU Symp. 113: Dynamics of Star Clusters, 

541-577 

Wilkinson, M. I., & et al. 2003, MNRAS, 343, 1025 



Woitas, J., Leinert, C, & Kohlcr, R. 2001, A&A, 376, 982 



